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TRANSFER FUNCTION ANALYSIS OF A FLUID CONTAINED IN A 
LONGITUDINALLY EXCITED THIN-WALL CYLINDRICAL TANK 

by Carl F. Lorenzo 
Lewis Research Center 

SUMMARY 

The transfer function for a constant- thickness, thin- wall cylindrical tank containing 
a fluid is derived. The fluid in the tank is assumed to have a constant ullage gas pres- 
sure and a pressure distribution initially linear with longitudinal position. The treatment 
is one-dimensional, and the effects of damping are considered. The transfer function de- 
veloped is the response of pressure at any longitudinal position to velocity at the tank 
base. General expressions for the resonant and nodal frequencies are given, and the 
time response to a sinusoidal driving velocity for an undamped liquid is derived. A nu- 
merical example which shows the general form of the transfer function solution is given, 
and a typical liquid- oxygen tank is analyzed. 


INTRODUCTION 

The occurrence of longitudinal oscillations (known as pogo) of significant amplitude 
in launch vehicles has produced efforts to define the problem analytically. This dynamic 
instability, which is usually due to the coupling of the vehicle structure (in a longitudinal 
mode) and the vehicle propulsion system, is known as the longitudinal structural dynam- 
ics problem. Analytical definition of the problem requires dynamic descriptions of the 
vehicle structure, the tankage characteristics, the propellant lines, the pumps, and the 
engines and their intercoupling. 

Studies of the overall problem (refs. 1 to 4) have placed emphasis on the modeling 
of pump- inlet and feed-line characteristics and other important areas, but have assumed 
relatively simple tank models for predicting pressure at the feed-line inlets. Generally, 
these studies have represented the propellant-tank characteristics as an equivalent single 
lumped mass or some minor modification of this concept. This study was conducted to 



determine the transfer function of the tank with the fluid mass considered to be distributed. 
Other studies (e. g. , refs. 5 and 6) have determined the modal characteristics of the pro- 
pellant tank. 

In this investigation, the propellant tank was isolated as a single component in the 
longitudinal system, and transfer functions for the tank (i. e. , response of pressure to 
driving point velocity) were derived in order to be useful either in a computer simulation 
of the entire system or in a direct calculation. 

The basic model considers the longitudinal wave propagation effects of a compress- 
ible liquid in a thin-wall, constant-diameter tank. For this model, the transfer function 
for the pressure at each longitudinal position to the velocity at the tank base is derived. 

The model used is one -dimensional and includes the effect of damping. In general, the 
problem is analogous to a modified form of the shorted load transmission line. 


SYMBOLS 


A area, sq ft 

a dummy variable defined by eq. (41) 

B fluid compressibility, sq ft/lb 
b dummy variable defined by eq. (42) 
c reciprocal of effective wave velocity, sec/ft 
D amplitude of sinusoidal input velocity 
d index for inverse transform, [n - (1/2)]tt/o! 

E modulus of elasticity, lb/sq ft 

f frequency, cps 

G dummy variable 

2 

g gravitational constant, ft/sec 
H dummy variable 

2 

K proportionality constant per unit length for fluid damping, lb-sec/ft 
K* (K/A dx)[B + (2R/Er)], sec/sq ft 
k constant pressure above liquid, lb/sq ft 
<£ Laplace transform operator 
S. height of fluid in tank, ft 
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M 


constant, lb/cu ft 
m, n index numbers 

P pressure, lb/sq ft 

R radius, ft 

s Laplace variable 

t time, sec 

u dummy variable for hyperbolic indentities 

V velocity, ft/sec 

y volume, cu ft 

v dummy variable for hyperbolic identities 

W weight, lb 

x distance coordinate, ft 

y dummy variable for inverse transform computation 

a constant 

e strain 

0(x,t) P(x, t) - P(x, 0), lb/sq ft 
X phase angle, deg 

p fluid weight density, lb/cu ft 

cr stress, lb/sq ft 

t wall thickness, ft 

<p(x,s) jar[e(x, t)] 
u> angular frequency, rad/sec 

Subscripts: 

FC fluid compression 

N nodal 

OV original volume 

R resonant 

S stored 

U undamped 

3 


AV wall volume change 
Superscript: 

differentiations with respect to time 
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(a) Forces. (bl Flows. 

Figure 1. - Coordinate system. 


DERIVATION OF BASIC EQUATIONS 

If the one-dimensional longitudinal motion of 
liquid in a cylindrical tank is considered, the co- 
ordinate system shown in figure 1 can be estab- 
lished. 

Considering the conservation of mass for a 
cylindrical element with area A and height dx 
gives 


8Wg 

at 





(i) 


The weight stored Wg in the elementary volume at time t is the sum of the 1 weight 
in the original volume, the weight stored in the change of volume of the walls, and the 
weight stored by fluid compression; hence, 


W S " W S, OV + W S, A V + W S, FC 


( 2 ) 


where 


W S, ov = pA * (3) 

and 

W S, FC = pABP * (4) 

The weight stored due to a small volume change is 

W S,AV = pAV (5) 

For a circular element of radius R, 

AV = 2ttR dx AR 
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If stresses remain within the elastic limit, the tank walls will obey Hooke’s law, or 


e = <Z_ _ AR 
~ E ~ R 

The hoop stress in the wall is a = PR/r; then, if the induced strain due to longitudinal 
stress (through the Poisson's ratio effect) is neglected, the radius change is given by 


AR = 


PR 

Et 


Hence, the stored weight due to volume change is 




( 6 ) 


Therefore, 


w g = + BP + ljpA dx 


(7) 


For constant-thickness walls and a constant-radius cylinder, this result can be substi- 
tuted into equation (1) to give 


\ et/ at ax 


( 8 ) 


If, instead, velocities are considered, 


\ Et / at ax 


(9) 


In order to derive the equation of motion, Newton's law is applied to an element acted 
upon by forces PA and - [p + (9P/9x)dx]A and a damping force. The tank base velocity 
V(0, t) is used as an approximation to the wall velocity (assumes longitudinal tank rigidity), 
and damping is assumed to be a linear function of the relative velocity. Further, since 
the tank wall is "thin” (t/R « 1), the associated mass can be neglected yielding the 
following equation of motion 
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( 10 ) 


PA - (p +i£ dx)A - W - k[v (x, t) - V(0, t)] dx = — — 


W dV 


dx 


g dt 


or 

[V(x,t) - V(0,t)] =0. (11) 

dx [a ) g at 

Differentiating equation (9) with respect to time and equation (11) with respect to x and 
eliminating the cross derivative terms gives 



c 2S?P + k.3P 


ar 


at 


( 12 ) 


where 



Figure 2. - Typical tank under consideration. 


and 


c 2 = p/b + 2*) 

g \ Et / 


> 





(13) 


When K* = 0, equation (12) is the undamped wave 
equation, and c is the reciprocal of the wave velocity. 
The term K* 3P/3t accounts for fluid dissipation losses. 


BOUNDARY CONDITIONS 

The system of interest and the boundary conditions 
are illustrated in figure 2. The ullage gas pressure above 
the liquid is assumed to remain constant. (The effect of 
acoustic variations in the ullage gas is given in ref. 7. ) 
Thus, 

P(£,t)=k (14) 

The response of the tank due to the sinusoidal motion of 
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a rigid tank bottom is desired; therefore, 


V(0, t) = D sin cot 

(15) 

Equation (15) is transformed into a pressure relation by considering equation (11) at 
x = 0: 

3P(o,t) ip _ p 3V(0, t) 
0X g 3t 

or 

(16) 

W(0.t) = . p .pDa, cosa , t 
dx g 

If a linear initial pressure distribution is assumed, 

(17) 

P(x, 0) = k -M(x - it) 

With the fluid at rest initially, 

(18) 

3P(x,0) = q 

at 

(19) 

Under normal conditions, M = p; however, in certain special cases ( 

[i. e. , bladder pres- 

surization systems for tanks under zero-gravity conditions (m = 0)), 
is desirable. The transformation 0 

the added generality 

0(x, t) = P(x, t) - P(x, 0) 

(20) 

is introduced into equations (12), (14), (17), (18), and (19), respectively, to simplify the 
initial conditions; then, 

d 2 e _ 2 d 2 e + K * de 
ax 2 " at 2 9t 

(21) 

0(£,t) = 0 

(22) 

cos wt 

dx g 

(23) 
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0(x, 0) = 0 


(24) 


jg (*,0) = Q 

at 


( 25 ) 


TRANSFER FUNCTION SOLUTION 


The Laplace transformation is applied to the preceding equations and boundary condi- 
tions. The transform of 0(x, t) is defined by 


<p(x, s) = jS?[0(x, t)] = J 0(x, t)e" st dt 


(26) 


where s is the Laplace operator. Since the initial conditions are zero (eqs. (24) and 
(25)), the transform of equation (21) is 


8 y( x ? s ) - (c 2 s 2 + K*s)<p(x, s) = 0 

,2 

ax 


(27) 


Equation (22) becomes 


(p{t,s) = 0 


(28) 


and equation (23) becomes 

d<p (0, s) _ M - p _ pDto s 

dx s g 2 2 

B s + a> 


(29) 


Since s is a parameter, equations (27) and (29) can be considered as ordinary differen- 
tial equations. 

The solution of the following equations and conditions is thus required: 


d A-gl - ( c 2 s 2 + K*s)(p(x, s) = 0 
dx 2 


(30) 


<p(a,s) = o 


(28) 


dcp(0, s) _ M - p _ pDo) 
dx s g 


s 



( 31 ) 
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The solution of the ordinary differential equation (30) can be written in the form 


cp(x, s) = G sinh ^Vc 2 s 2 + K*s + H cosh ^Vc 2 s 2 + K*s x^ 
Then, applying condition (28) yields 


G sinh (V c 2 s 2 + K*s f) + H cosh(Vc 2 s 2 + K*s = 0 
and applying equation (31) gives 


d< ?(°> s) e G |/ C 2 S 2 + k »b = M - p _ pDo) s 

dx s e 2 2 

» S + 0) 


from which 


G = 


M - p pDw 


J 2 2 

ic s + 


K*s s 


l/ 2 2 
T C S . + 


K*s (s 2 + w 2 ) 


Equations (33) and (35) combine to give 

H = -G tanh (Vc 2 s 2 + K*s i) 


( 32 ) 


(33) 


(34) 


(35) 


(36) 


The solution of equation (30) is thus equation (32), where G and H are given by equa- 
tions (35) and (36), respectively. After some manipulation, equation (32) can be written 
as follows: 


<p(x,s) = - 


M - p pDfa) 


J 2 2 

Tc s + 


K*s 


g 


sinh 


V 2 2 

f c s + 


K*s (4 


Vc 2 s 2 + K*s (s 2 + w 2 )J cosh (Vc 2 s 2 + K*s 



(37) 

In order to evaluate the transfer function, that is, the ratio of the Laplace transform 

of the system response to the Laplace transform of the driving function with quiescent 

2 2 

initial conditions, M is set equal to p. Since V(0, s) = Du/(s + w ), 
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s 


( 38 ) 


<?(x,s) _p 

V(0, s) g 


s fsinh [Vc 2 s 2 + K*s (SL - x)J 
Vc 2 s 2 + K*s cosh (Vc 2 s 2 + K*s ij 


This function is the transfer function sought in terms of the Laplace variable s and the 
system parameters. The form is such that it can be substituted into many existing simu- 
lations that derive tank pressure to driving point velocity by using a single lump approxi- 
mation. 

In some situations, having the time response of this system is desirable. Hence, 
the time response is determined in the appendix for the case with zero damping and gen- 
eral linear initial conditions ompressure distribution. The pressure response in time 
due to a sinusoidal input in tank base velocity is given by 


0(x,t) = iT 1 [^(x, s)] = -(M - p) 


00 

(£ - x) + — — \ izll— sin 6c(SL - x)cos 6t 

2„ Z-f *2 

c 1 fcl 6 


/pDco\ _ x ) cos \ ' ill) 1 S in[6c(£ - x)](6 2 cos 6t - co 2 cos cot) ► 

v */] ch ^ 6 2 <o 2 - 6 2 j 


where 
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JL 

cSL 


(A8) 


(A9) 


The amplitude ratio and phase relations corresponding to equation (38) are important 
for experimental testing and for various other analyses. Substituting ico for s to ob- 
tain the magnitude and phase relations associated with steady- state vibration gives 


<?(*, M _ P 

V(0, ico) g 


ICO 


sinh 


J 2 2 

f -c co + 


K*coi (SL - x) 


V-c 2 co 2 + K*coi 1 cosh ^V-c 2 co 2 + K*coi SL ) 


(39) 
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m 1 

The quantity \-c + K*ct>i can be considered a complex number; hence, if 
V-c 2 w 2 + K*wi is replaced by a + ib, equation (38) can be written 


where 


<?(x» io)) _p j 

V(0, io>) g a + ib 


sinh[(a + ib)(£ - x)] 


cosh [(a + ib)f] J 


( 40 ) 


and 


a 


= _Va> 


4 + (K*co) 2 - cV 


b = y i [Vc 4 co 4 + (K*w) 2 ' + cV 


(41) 


(42) 


Equation (40) and its parameters a and b are similar in form to the shorted load trans- 
mission line equations. They differ basically in the form of dissipation considered; that 
is, for transmission line equations, damping depends on absolute velocity (current) in- 
stead of the relative velocity considered here (i. e. , eq. (10)). 

The amplitude ratio can be obtained from equation (42) as 


(43) 


<?(x, iw) 


pcoi 

|sinh[(a + ib) (£ - x)] | 

V(0, io>) 


g(a + ib) 

|cosh[(a + ib)£]| 


The following relations are substituted into equation (43) : 


pan 


g(a + ib) 


_ pw 
g 




| sinh(u + iv) | = ^jsinh^u + sin 2 v 
| cosh(u + iv) | = Vsinh 2 u + cos 2 v 


i 2 + b 2 = Vc 4 o> 4 + (K*gl>) 2 ' 


(44) 

(45) 

(46) 

(47) 
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which yield as the expression for the amplitude ratio 


<p(x, io>) 


po) 

-t i 

V(0, iw) 

g 

[c 4 u> 4 + (K*w) 2 ] 

1/4 y 


sinh 2 [a(C - x)J + sin 2 [b(£ - x)] 
sinh 2 (al) + cos 2 (bf) 


(48) 


Thus, since the amplitude ratio has been determined, the phase relation is determined 
from equation (40) and from the relations 


sinh(u + iv) = sinh u cos v + i cosh u sin v 
cosh(u + iv) = cosh u cos v + i sinh u sin v 


(49) 

(50) 


Then 

<p(x, io>) _ pun / sinh[a(f - x)]cos[b(£ - x)] + i{cosh[a(£ - x)]sin[b(£ - x)]|' 

V(0, iw) g(a + ib) y cosh(a£)cos(M) + i[sinh(al)sin(bf )] 


(51) 


Since arg(z j/zg) = arg z^ - arg Zg, the general phase angle solution X (in deg) is given 
by 

X = 1§2 /tan _1 |coth[a(£ - x)]cot[b(£ - x)]| - tan - 1 [tanh(af )tan(bl)] - tan"*^:j\ (52) 


For the special case of zero damping (K* = 0), the amplitude ratio (eq. (48)) reduces to 


(53) 


(p(x, iw) 


sin[ct>c(£ - x)J 

V(0, ico) 

u gC 

cos(wcf) 


or 


<p(x, M 


|V(0,iw) 


= -£- | tan(o> c£ ) cos (w cx) - sin(o)cx) 


( 54 ) 


U 


gc 


The phase angle for the undamped case is determined from equation (40) for K* = 0: 
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( 55 ) 


<p(x,io)) _p fsin^cd - x)]l . 
V(0, icu) gc\ cos(wc£) J 


From equation (55) it is clear that 


\j=< 


+90 


-90 L 

v. 


sin[a>c(£ - x)1 

for ^ > 0 

cos(cocf) 

sin[u>c(£ - x)] 


for 


cos(wcfi) 


<0 


(56) 


The driving point impedance relations can be determined in a similar maimer by 
letting x = 0 in equations (48) and (52). The resonant and nodal frequencies can be cal- 
culated for the undamped case by considering the expression for the amplitude ratio 
(eq. (54)). 


Resonant Frequencies 

For the undamped resonant frequencies, the amplitude ratio becomes infinite. If the 
bottom of the tank (x = 0) is considered, from equation (53) 

|tan(ci>c£) | = 00 


which requires 


and 


= 



(rad/sec) 


n = 1, 2, 3, • • • 



n = 1,2,3, * • • 


(57) 


(58) 


This result can also be shown by a consideration of the time response (eq. (A8)). 
This expression for the resonant frequency has the same form as that for a closed organ 
pipe, as might be expected, since the undamped case was considered. 

Equations (57) and (58) are not strictly applicable for all tank stations since certain 
x values produce the indeterminate form zero over zero, which when evaluated can be- 
come a finite amplitude ratio. This fact is, however, of no practical consequence. 
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Nodal Characteristics 


The location of the nodes (the condition under which the amplitude ratio becomes 
zero) for the undamped case, is determined from equation (53), hence 

sin u>c(£ - x) = 0 

which requires 

x = L - — for 0 £ x s £ n = 0, 1, 2, • • ■ 

OJC 

For a given x -position the frequency at which nodes will occur is given by 


w, T = — — for 0^x<i n = 0, 1, 2, • • • 

w c(£ - x) 

At resonance, the nodes will be determined by the formula 


n 


x 



m 


l 



for 0 < x ^ £ n = 0, 1, 2, • • • m = 0, 1, 2, • • • , n - 1 


(59) 


(60) 


NUMERICAL EXAMPLE 


Because of the relative complexity of forms of the transfer function (amplitude ratio 
and phase), a special case was evaluated to indicate the form of the solution. For sim- 
plicity, the following tank and fluid characteristics were assumed: SL = 1. 0, p/g = 1. 0, 
and c = 1.0. Under these conditions, the magnitude and the phase angle were evaluated 
from equations (48) and (52), respectively. 

For the undamped case (K* = 0), the effect of the x position in the tank is given in 
figure 3. The bottom of the tank (x = 0) shows nodes at co = ror (for n = 0, 1, 2, • • • ), 
while for x = 0.25 and x = 0.75 the nodes appear at 4njr/3 and 4mr, respectively. 

The amplitude ratio plot has a period of ir radians per second for the x = 0 position 
and a period of 4ir radians per second for x = 0. 25 and x = 0.75. This periodicity of 
the magnitude with co (or, in general with cod) occurs when an integral number of half- 
cycles of the numerator and an integral number of half-cycles of the denominator (of 
eq. (53)) occur at a common value of coc£ . For x = 1/3, for example, 


(p(x,ia>) 

= _P 

sin^ coc&J 

V(0, ico) 

gc 

cos(cod) 


( 61 ) 
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Figure 3. - Effect of position on response of pressure to velocity for undamped case 
(magnitude). 


When <ji c£ = 3ir, two half- cycles of the numerator and three half- cycles of the denomina- 
tor will have occurred and this sequence will be repeated for each successive 37r inter- 
val of wc S.. In general, the amplitude ratio for zero damping is periodic with o>c£ 
and has a period of ran, where m is the denominator of the least common fraction 
(£ - x)/£ for a given x position. 

The resonant frequencies for all the locations shown occur at [n - (1/2 )]jt (for 
n = 1, 2, 3, • • •)> as predicted by equation (57); however, this is not the case for all 
locations, as discussed in the section Resonant Frequencies. Since the undamped case 
was considered, the phase angle is either 90° or -90°, as determined by equation (56), 
and is therefore not plotted. 

The effect of damping was studied, and typical results are shown in figure 4. The 
position x = 0 (line inlet) was chosen because of its special importance in many cases. 
The amplitude ratio plot (fig. 4(a)) indicates that increasing damping K* decreases the 
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Phase angle, \ deg Amplitude ratio, |<p(x, M/V(iu)|, Ib-sec/cu ft 



Figure 4. - Effect of damping on response of pressure to velocity 
at tank base. 


height of the resonant peaks as would be 
expected. The shifting of the resonant 
frequency, however, appears to be negli- 
gible. The effect of damping on phase 
angle (fig. 4(b)) also shows negligible 
shifting of the zero crossing at the reso- 
nant frequencies for the cases considered. 
The phase angle at zero frequency ap- 
proaches 90° for the x = 0 position in- 
stead of zero that might be expected, be- 
cause of the form of the damping assumed 
in equation (10). 

LIQUID-OXYGEN TANK ANALYSIS 

In order to interpret the results of 
this study in a quantitative manner, the 
dynamics of a tank having missile dimen- 
sions is considered. The assumed tank 
has a 5-foot radius and a 0. 0173-inch 
average wall thickness and is filled to a 
depth of 20 feet with liquid oxygen. For 
stainless -steel walls, these dimensions 
imply (eq. (13)) an effective acoustic ve- 
locity 1/c of 487 feet per second, or a 
value of c of 0. 00205 second per foot. 
The amplitude of the resonant peaks is, 
of course, dependent on the damping mag- 
nitude. Since experimental values of K 
or K are unknown, a plot of the effect 
of K* on peak amplitude is presented in 
figure 5. Over the range shown, peak 
amplitude is nearly linear with K* on the 
log- log plot with a slope approximately 
equal to 1. For an assumed value of 
K* = 2x10" the model indicates that the 
amplitude ratio | <p/V | at the first reso- 
nant peak is 1. 14x10^ at the tank base. 
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mi 

'I 

?• 



Figure 5. - Effect of damping on amplitude ratio of first resonant peak of liquid-oxygen 
tank. Effective wave velocity, 487 feet per second; fluid weight density, 71. 14 pounds 
per cubic foot; height of fluid in tank, 20 feet 


Hence, for an input velocity amplitude of 0. 2 feet per second, the corresponding pressure 
amplitude is 15. 8 pounds per square inch. Pressure amplitudes of this magnitude have 
been observed in unpublished data. 


SUMMARY OF RESULTS 

The transfer function of a fluid contained in a longitudinally excited tank was devel- 
oped, based on the following assumptions: 

(a) Constant ullage pressure 

(b) Constant-thickness, thin, cylindrical walls 

(c) Negligible Poisson's ratio effect in the tank walls 

(d) Initial pressure distribution linear with position 

(e) One-dimensional motion (no geysering or three dimensional effects) 

(f) Linear dissipative losses 

The following determinations resulted from the analysis: 
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1. The transfer function (pressure to tank base velocity) was derived in terms of the 
Laplace variable and in terms of the amplitude ratio and the phase angle. The amplitude 
ratio and phase angle for undamped motion were also determined. The special case of 
driving point impedance is readily determined in the preceding cases by allowing x = 0, 
where x is the distance from the tank base. 

2. The time response in pressure to a sinusoidal input velocity at the tank base was 
determined for the undamped case. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 9, 1966. 
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APPENDIX -TIME RESPONSE 


The time response of the system under study is of some interest in understanding the 
physics of the problem. The undamped response of the system is considered. In this 
case, from equation (37), the inverse Laplace transform of the following equation is re- 
quired: 


<p(x,s) = - 


M - p _ pDio 
.2 gc 


cs 


1 \ fsinh[cs(f - x)]l 

+ w 2 /t cosh ( cs£ ) J 


(Al) 


By reference 8, 


g>-l( sinh ys 
Vs 2 cosh as) 




il 

a 


for 0 < y < a 


(A2) 


also, 


<T X 




= — sin o>t 


OJ 


(A3) 


The convolution theorem (ref. 9) is applied to these functions to obtain 



Performing these integrations and manipulating trigonometric forms produces 
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<e 


-l 


sinh ys 


2 2 2 / 
s + co s cosh as/ c o' 


— (1 - cos cot) 
2 


where 




sin yd(cos td + cos 


cot) 


(A 5) 



(A6) 


2 

Suppressing the s term in the denominator by differentiating the time result twice with 
respect to t (ref. 10) determines the inverse transform: 

Jtahyi) = y cos cot 

V+CO 2 c0sh °7 

2 2 

sin yd (d cos td - co cos cot) (A7) 



By application of the inverse transforms (eqs. (A2) and (A7)) to equation (Al), the time 
solution is determined as 


0(x,t) = ir^cpix, s)] = -(M-p) 


(1 - x) + • 


y (ii^si 

,2 


sin 6c(£ - x)cos 5t 


cl w 6 


+ J(£ - x)cos cot — — sin[6c(£ - x)](6 2 cos 6t - co 2 cos cot)! 

' ® / I c 2 £ ^—4 5 2 co 2 - S 2 I 

n=l 


(A8) 


where 
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(A9) 



l 

I' 
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